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I. INTRODUCTION 


In modern missile warfare, new technologies are producing faster, more 
accurate anti-ship cruise missiles (ASCM). They are able to fly at very low 
altitudes and perform a high speed ‘pop-up’ maneuver just before impact, in 
order to defeat the ship's missile defenses, as shown in Figure I. In order to 


counter the advanced ASCM's, new control schemes must be developed to 





Figure 1. Anti-ship cruise missile terminal maneuver 


enable the ship's surface-to-air missiles to counter this threat. 

The aim of this thesis is to illustrate the concept of Bang-Bang control and 
how it applies to control of vertically-launched surface-to-air missiles. Bang- 
Bang control can be used in the missile boost phase in order to turn the missile 
downrange in minimum time, and also in the terminal phase, where the missile 
must react very quickly in order to intercept a maneuvering target. 

This work applies the minimum-time (Bang-Bang) control concept to 


vertically-launched surface-to-air missiles. In Chapter II, the concept of 


control is defined and solved for a second order system. In Chapter III, the 
equations of motion for a vertically-launched missile in the boost and terminal 
phases are developed. In Chapter IV, the boost simulations comparing missile 
trajectory versus maximum thrust vector control (TVC) angle are presented. 
In Chapter V, the terminal phase missile simulations are presented for a non- 


maneuvering target and also for a maneuvering target. 


Il. CONCEPT OF BANG-BANG CONTROL 


In this section the form of the optimal control for a particular class of 


systems will be determined by using Pontryagin's Minimum Principle. 


A. PONTRYAGIN'S MINIMUM PRINCIPLE 


It shall be assumed that the state equations of the system are of the form 


x= Ax+Bu, | =i) 


where A is an n by n array, B is an n by m array, and A and B may be 
explicitly dependent on the states and time. For this research it will be assumed 
that there is a single input; therefore, m is 1. It is specified that the admissible 


controls must satisfy the inequality constraints 
N.s u(t) Ss N, (2-2) 


where N. and N, are the known lower and upper bounds for the control input. 
It is desired to drive x from an initial state x(ty) to to a desired final state x(t), 

where tg is the problem start time and t; is the problem end. The optimal control 
is that control which drives the state from its initial state to desired final state 
using the least "cost". The cost can be any desired measure of system 
performance. In the case of a minimum-time controller, the cost function can 


be represented as 


tf 
J= | dt=t,-ty. 
« (2-3) 


Pontryagin's Minimum Principle states that the optimal control, u*, which 


minimizes the cost function, must minimize the Hamiltonian, which is defined as 
H(x(t), u(t), p(t), t= 1+p™(t)(Ax(t)+Bu(t)) (2-4) 


where p(t) is a Lagrange multiplier vector, and is arbitrary. p(t) can be defined 


as the 'costates' of the system, and p(t) can be written as 


oO 


2-5 
P2(t) — 
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The necessary condition for the optimal control u* to minimize the cost 


function J is 
H(x*(t),u*(0),p*(t),t) S$ H(x*(0),u(t),p*(t),t) (2-6) 


for all times between to and tsand for all admissible controls. The asterisks in 


equation (2-6) represent the optimum values, and thus, 
1+p*T(t)Ax*(t)+p*!()Bu*(t) 21+p*T(t)Ax(*t)+p*!(t)Bu(t) (2-7) 
Therefore 
p**(t)Bu*(t) Ss p*1(t)Bu(t) (2-8) 


for all admissible u(t) and for all times from t.,;;, to tsa, If u(t) 1s constrained to 


+N, the optimal control U*(t) is 


* —N f : >0 * 
lS e P2 4 -Nsign(p; | (2-9) 
+N for p, <0 


This indicates that the time-optimal control is ‘Bang-Bang’; that is, the 
optimum control switches between its maximum positive and negative values. 
There are also three ideas which deserve mention: 

(1) (Existence) 
If all the eigenvalues of A have nonpositive real parts, then an optimal 
control exists that transfers any initial state xo to the origin. If there 
are positive eigenvalues (i.e. unstable roots), there may be some 
region of the state space where the system is uncontrollable. 
(2) (Uniqueness) 
If an optimal control exists, then it is unique. 
(3) (Number of Switchings) 
If the eigenvalues of A are all real, and a (unique) time-optimal 
control exists, then the control can switch at most (n-1) times. 
Thus, an nth-order system having all real, non-positive eigenvalues has a unique 
time-optimal control that switches at most n-1 times (not counting switching off 
at time trina). For complex conjugate eigenvalues, more than n-1 switchings 


may be required.[1] 


B. APPLICATION TO A SECOND ORDER SYSTEM 
In order to illustrate the minimum principle, consider the second order 


differential equation 


Cae (2-10) 


or in state space, 


ie +] ju. 
©) QO O X95 1 (2 ] 1) 
From the minimum principle, the time-optimal control for this system is £N. 


Thus the segment of optimal trajectories can be found by integrating the state 


equations (2-11), with u=4N, from time tp to t. This yields the results 


X,(t)=ENt+C, (2-12) 
2 
AE) ee cee (2-13) 
where 
GC = +Ntp + X(t) (2-14) 
tea 
C, = FN —-_-C,t, +x, (t, ). 
and : De la aa (2-15) 


C, and C, are functions of the initial conditions and can be treated as constants. 
Time can be eliminated from equations (2-12) and (2-13) by squaring the first 
equation, dividing the result by 2N, and comparing the result with equation 


(2-13) to obtain 


n (= x(t) + Cy for U=+N (2-16) 


x(t) =x (D+Cy for U=-N (2-19) 


where 








2N (2-18) 
2N (2-19) 


C; and C, are functions of the initial conditions and may also be treated as 
constants. Equations (2-18) and (2-19) each define a family of parabolas, as 
shown in Figures 2 and 3. 

Analysis of Figures 2 and 3 reveal the controls corresponding to the 
following situations: 

1. u*=+1 implies that the initial state xo lies on segment AO at time to. 

2. u*=-1 implies that the initial state xo lies on segment BO at time to. 

This defines the optimal control to reach the origin. In order to reach these 


zero trajectories, the opposite value of control is required. It can be seen from 


Figure 2. Trajectories for u=-N 


= 


Figure 3. Trajectories for u=+N 


| > 


Figure 4 that the two segments AO and BO form a single curve. This curve is 
known as the switching curve. All states that fall above the switching curves 
will result in the control input U=-N, and those states which fall below the 
switching curves will result in the control input U=+N. This curve can be 


represented as 


X,(t) = - x, (tx (0) (2-20) 


Thus, the control law can be obtained by moving all the terms of equation 


(2-20) to one side: 


U(t)=-Nsign{ «,(0) +5 x9(0fx2 (0) (2-21) 


ee — = 





Figure 4. Optimal switching curves 


Examples of the use of this switching law is shown in Figure 5. From any 
arbitrary initial state, the trajectory will follow the parabola that passes through 
the start point until the trajectory reaches the optimal switching curve. Then the 
control switches signs, and follows the optimum parabolic path to the origin. 

Equation (2-21) is the mathematical representation of the Bang-Bang 
control law. This control law can be used in numerous applications where it is 
desired to transfer a state from initial conditions to desired final conditions in 
minimum time. Two examples that are of interest in this research are: 

1) drive a state, such as a position, to zero in minimum time. The control 


law for this portion of the missile flight is of the form 


UW) =-Nsiga{ x (1) + x2(x2(0)} (2-22) 


oN 
Na 


Figure 5. Optimal control switching schemes 


where x,(t) is the position at time t, x(t) is the velocity at time t, and N is the 
maximum admissible force that can be used to accelerate or decelerate the 
state; and 
2) drive the time rate of change of a state to zero. In this case it is 
desired to keep a state constant, as opposed to driving it to zero. The form of 
the controller is also given by equation (2-22), but x, in this case will represent 
velocity and x2 will represent acceleration, with N being the maximum 
admissible control with which to change the states. 
In the next chapter, the application of the second order Bang-Bang 


controller to boost phase and terminal phase missile control will be examined. 





Hl. DEVELOPMENT OF EQUATIONS OF MOTION FOR A 
VERTICALLY-LAUNCHED MISSILE 


The equations of motion for a missile can be derived from Newton's second 
law of motion, which states that the summation of all external forces acting on a 
body must be equal to the time rate of change of its momentum, and the 
summation of the external moments acting on a body must be equal to the time 
rate of change of its moment of momentum (angular momentum). The time 


rates of change can be expressed by two vector equations: 


sn= Sa) on 
and 
OM, = (3-2) 


where 1 indicates i'4 cartesian coordinate of the vector with respect to inertial 
space. By definition, H is the angular momentum, or moment of momentum, of 
a revolving body.[2] 
There are several assumptions that are made in order to simplify the 

problem: 

1. The mass of the missile remains constant. 

2. The missile airframe is a rigid body. 

3. The earth is an inertial reference, and unless otherwise stated, the 


atmosphere is fixed with respect to earth. 
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In order to completely describe the motion of a missile in three dimensions, 
a total of six nonlinear differential equations must be solved. However, certain 
assumptions will be made in order to reduce the number of equations needed to 
adequately describe the motion of the missile: 

1. The X axis will be assumed to be downrange from the launch 
platform to the target, the Y axis will be crossrange, and the Z axis will be 
altitude. 

2. The type of missile studied here is known as a skid-to-turn missile, 
as this missile uses direct side force to turn. Thus, there will be no coupling 
terms between the pitch and yaw axes, and one set of equations will describe 
the motion in each of the Y and Z planes. 

3. It will be assumed that the missile is roll stabilized, 1.e., that there 
will be no turning moments about the X axis. 

4. Since there are no cross-coupling terms between the pitch and yaw 
axes, the simulation can be broken into two problems, one where the pitch 
angle is held constant, and the other where the yaw angle is held constant. In 
this discussion the yaw angle will be assumed to be constant at zero. Using 
these assumptions, a 3 degree-of-freedom model can be developed. The force 
diagram of a 3 DOF missile in flight is shown in Figure 6. The governing 


equations of motion of this missile are: 


DEE =mxk= -Fi Siny, + Tcos(0,, +9) - Fig COS Ym (3-3) 
> FE, =mz=-—mg + Fig cosy, + Tsin(0,,+5) - Fuse siny, (3-4) 


De Meg = Ieg9 = — Fig pg COSY m — Fane pe Si eet (3-5) 
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Figure 6. Forces acting on a missile in flight 


where: 
m missile mass (Kg) 
g gravity (9.81 m/s2) 
Fiirt aerodynamic lift force (N) 
Farag aerodynamic drag force (N) 
Ym angle of velocity vector from reference 
8., angle of missile head from reference 
6 TVC deflection angle 


I-g moment of inertia of the missile about the center of gravity 
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Lpg distance from aerodynamic center of pressure to center of 
gravity 
Lig distance from TVC control to center of gravity 
The aerodynamic force acting on a lifting body, in this case the wings and 
control surfaces of the missile, may be resolved into two components, which 


are 


Fig =Cy 2 v2s (3-6) 


and 


P x72 
F,.2 =Cp~V,9 
te (3-7) 
where (p/2)V,,2 is the dynamic pressure in N/m2. 
The lift and drag coefficients are functions of angle of attack (a), and Mach 
number. Actual values for C, and Cp must be measured, as they are highly 
dependent on the airfoil geometry. However, Cy, and Cp may be modeled 


satisfactorily using the following approximations: 





C, =0.la (3-8) 
— 2T 2 
Cy = D _S (1 +0.20 i) (3-9) 


These coefficients are in reasonable agreement with reference [2] at low 
angles of attack (a<10°). The only time the missile can expect to exceed 10° 


angle of attack is during the tip-over phase, when velocity, and therefore lift and 
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drag forces, are small compared to the thrust force. This is considered an 


appropriate approximation for the scope of the model under discussion. 


A. BOOST PHASE CONTROLLER 

The primary objective of the missile control system in the boost phase of 
flight is to tip the missile over and point the missile in the general direction of the 
target in minimum time. Therefore, 9, the angle that the missile makes with the 
horizontal reference, is a logical choice for a control variable. If the force that 
creates the moment that tips the missile over can be controlled, a general 
control algorithm can be formulated. This is shown graphically in Figure 6. 

At launch, the missile is assumed to be in a vertical position, with the 
velocity vector pointing straight up, i.e. 8 = 90° and y = 90°. In the initial boost 
phase, it is assumed that the first two terms of equation (3-5) are negligible as 


compared to the third. Therefore, equation (3-5) may be simplified to 
- —TL,,sind 
9 = —_—_>— 
I (3-10) 


cg 


In state space, this can be represented as 
QO 116} |O} -TL,, si 
Pie 7}4{0 | —thesine (3-11) 
a QO 0] 8 l a 


The only term in equation (3-11) that can be used to control @ is 6, the TVC 
deflection angle. A simple and effective method for controlling missile tip over 
is to employ the Bang-Bang controller, developed in Chapter II, to drive 8 to 0 
(horizontal flight) in minimum time, as shown in Figure 7. Thus, the control law 


iS 


VS 
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Figure 7. Boost Phase Diagram 
2 1 so. (3-12) 
u(t) =-Nsign{ 0(1) + 6(0]6¢0) 


where N is the maximum torque available to rotate the missile and is 


determined from the relation 





de 
N==——= sin(d 
cg 


ae (3-13) 
and Omax is the maximum thrust vector control (TVC) angle available. 


This indicates that for minimum tip-over time, dnax=90°, However, this 


may not be desirable since all of the thrust will go to rotate the missile, while 
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none of the thrust will accelerate the missile downrange. This amounts to the 
missile 'pin-wheeling' in the air. Therefore, TVC angles of less then 90° should 
be used. Several TVC angles will be examined to show the tradeoffs involved. 

When the missile achieves horizontal orientation (i.e. 82=0°), the control law 
will continue to switch from +N to -N in an attempt to keep the missile at 
exactly 0°. This phenomenon is known as ‘chattering’. At this point it is no 
longer desirable to use maximum control effort to correct for small changes in 
8. A satisfactory method for removing the chatter from a Bang-Bang 
controller is to employ linear control in a small region, of width €, about the 
desired value, as shown in Figure 8. 


This adapted Bang-Bang controller can be expressed mathematically as 


PE 





ieee il eure 
Nsign( 0(t)4 a ajo] act) e 6(t)/6(t)| 
u(t) = 
-2(0c0) +5 6(1]6¢0) ~¢ < (t) + 6(t)/6(1)|<e (3-14) 
| ¢€ 2N 2N 





Figure 8. Bang-Bang Control Law with Linear Zone 
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B. TERMINAL PHASE 
1. Missile-Target Geometry 
In this section, it will be assumed that the vertically launched surface- 
to-air missile has completed tip-over, and has settled out to a constant speed 
and altitude. The start point for the terminal engagement of the target begins 
with the missile seeker acquiring the target. The geometry of the missile-target 


engagement is shown in Figure 9. 


Missile 





Figure 9. Missile-target geometry 


In describing the geometry of the missile-target encounter, several 
parameters must be calculated. The angle the missile's velocity vector makes 


with the inertial reference frame is defined as y: 


y = tan! Vym 
m V ; (3-15a) 
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The magnitude of the missile velocity vector V,, is given by the relation: 


V, = V2, + V2, (3-15b) 


where Vym 1s the vertical component of the missile velocity vector and Vxm is 
the horizontal component of the missile velocity vector, measured in inertial 
coordinates. 

Similarly, the magnitude of the target velocity vector and and angle it 


makes with the inertial reference is given by 


oe 2 z 
V= Vin -} Wye (3-16a) 


and 


val (3-16b) 


In computing the variables used to calculate the line of sight and its 
derivatives, it is convenient to describe the relative position, velocity, and 


acceleration between the target and the missile: 


XP Sn (3-17a) 
a= Yn . (3-17b) 
We xt - Voan (3-17c) 
= eave (3-17d) 
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As = Axt oe Ag (3-17e) 


The line of sight between missile and target is given by O, and is 


defined as: 


va ( 
X (3-18) 
where X and Y are the relative cartesian positions as defined by equations (3- 
17a) and (3-17b). The first and second time derivatives of 6, denoted G and 6, 


are computed using the relations 


__(XV,-YV,} 
on R? (3-19) 


and 


._XA,-YA, 2(XV,+¥V,)(XVy - YVx} 


G @ ad (3-20) 


In an actual missile-target engagement, it is very unlikely that the 
target line of sight rate and acceleration will be measured accurately enough in 
order to use G and go directly in the guidance algorithm. An alternative 
method is to calculate estimates of G and 6G, denoted by G and G , based on 
measured values of 0. This is accomplished through the use of a Luenberger 


Observer, as shown by the signal flow graph in Figure 10. 
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Figure 10. Luenberger Observer for estimating 6 and 6 


This observer can be described in state-space form by the relation 


0 


G6) |7ks 1 6] [k, 
A —2R a 
g — —k, jp” I 2 - k, O (3-21) 


where R is the range from the missile to the target, R is the range rate, and ky, 
k>, and k; are gains to be determined by the designer in order for the estimated 
States to closely follow the true states. 
2. Missile Guidance Laws 

The objective of tactical missile guidance is to keep the line of sight 
angle, 0, between missile and target constant. This is desirable because if the 
line of sight remains constant, while the range from missile to target decreases, 
an intercept will occur, as shown in Figure 11. Two missile guidance laws will 
be presented: classical proportional navigation and minimum time (Bang-Bang) 


control. 
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O represents line of sight angle 
between missile and target. 





Cf Missile 


Figure 11. Constant Line-of-Sight Intercept 


a. Proportional Navigation Control 

The method most often used in controlling homing missiles is 
known as proportional navigation, where the magnitude of the transverse 
control force is proportional to the rate of change of the line of sight (6). 
Lukenbill [3] conducted research into proportional navigation control and his 
derivation of the proportional navigation controller is given below. 

Figure 12 depicts the basic proportional navigation scheme. 
Assuming that the seeker head of the missile follows the target, the transverse 
acceleration perpendicular to the line of sight will equal the acceleration of the 


R vector in that direction. Mathematically, the acceleration of R is 


Ag =(R+@X@XR)ig +(2@XR+@XR)ig (3-22) 


where 


R ae missile/target line of sight vector 
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R = closing rate along R 
R = acceleration along R 
= angular rate of change of R in inertial space 
Ar = overall acceleration of R. 
At this point, a missile acceleration, A,,, equal to the target 
acceleration, A,, will make the line of sight parallel to its original direction. As 
long as R remains along R (@=0) a missile/target intercept is assured. So, the 


transverse acceleration command is 


A,-A, =@xR+2(@xR). (3-23) 


Ze 





Figure 12. Vectorial Proportional Navigation Scheme 


Assuming the line of sight rate is equal to the angular rate of change of R in 


inertial space, equation (3-23) now becomes 


A, = - RG+2RG6 (3-24) 


where 6 equal to @, and oO is equal to o. 

In the classical proportional navigation scheme, the missile course 
is one in which the rate of change of the missile heading is directly proportional 
to the rate of rotation of the line of sight vector from the missile to the target. 
As a result, this course change is intended to counteract the rotation of the line 


of sight, thus returning to a constant bearing course. The movement of the 
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missile and target cause the line of sight to rotate, resulting in a differential 
displacement between the missile and the target perpendicular to the range line. 
Figure 13 depicts this geometry. The proportional navigation guidance law 
attempts to generate an acceleration command, A,, perpendicular to the line of 
sight. 

Assume a gyro stabilized seeker head, as in the Sidewinder 
missile. If there is no torque applied to the gyro, the seeker will not rotate. 
Assuming the seeker tracks the target, the gyro angle will follow the line of 


sight. Applying the equation of motion for a gyro stabilized seeker 
=l@Q (3-25) 


where 
= applied torque 
= spin angular velocity 
I = moment of inertia of the gyroscope 


Q rate of precession of the gyroscope. 


Applying this to the case when the seeker head tracks the target, Q is then 


replaced by the rate at which the gyro is torqued in space. This is simply 0, 


Vertical Displacement 





Missile 


Figure 13. Missile Acceleration Orientation 
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Figure 14. Missile Acceleration Relationship 


which is the line of sight rate. Thus, equation (3-25) becomes 


L=Iqwo. (3-26) 


This torque is in turn applied to the control surface of the missile leading to the 


relationship 


A,, =kL=klo6 (3-27) 


where k is a constant of proportionality. Referring to Figure 14, a relationship 1s 
determined for A,, in terms of the rate of change of the missile flight path angle, 
Y,- Given the missile velocity vector at some point in time, V,,(t), and suppose 
the missile undergoes an acceleration, A,, during an interval of time, dt. The 
velocity vector is then displaced and is represented by the vector V,,(t+dt). The 
angle the vector is traversed is simply dy,,, the differential missile flight path 
angle. For small angles (which are guaranteed by making dt small) the 


following relationship is obvious: 


A_dt=Vidy_. (3-28) 
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Dividing equation (3-28) by the time interval, dt, the missile acceleration is 


defined as 


Ag = Va 2 = Vr 


Combining equations (3-27) and (3-29) 


V, 


m 


Ym = Kloo. 


Dividing through by Vy, the proportional navigation law becomes 


, _(Ko\, 
Ym V. 


Or 


Ym = NG. 


Ym, 


(3-29) 


(3-30) 


(3-31) 


(3-32) 


Equation (3-32) represents the classical proportional navigation equation where 


= rate of change of the missile heading 
Oo = rate of change of the line of sight 
N = proportional navigation ratio. 


The navigation ratio determines the sensitivity of the missile 


system. A high navigation ratio will lead to rather high gains resulting in large 


missile commands for small changes in the line of sight rate. On the other hand, 


small values for N will lead to small missile commands for a given 6. Larger 


navigation ratios are preferred for head on engagements and smaller ones are 


preferred for tail chase cases. For this research the navigation ratio is taken to 


be four. 


a 


In tactical radar homing missiles using proportional navigation 
guidance, the seeker provides an effective measurement of the line-of-sight 
rate, and a Doppler radar provides closing velocity information. [4] 

b. Bang-Bang Control Minimizing Line-of-Sight Rate 

In this control law, the goal is the same as in proportional 
navigation control: to drive the line-of-sight rate to zero, or in other words, to 
keep the line of sight constant. 


Mathematically, this control law can be stated as. 


u(t)=—-N sin( +51) (3-33) 


In order to alleviate chattering, a small linear region will be introduced, similar 


to the boost phase controller: 


ie . oe. . 
-Nsign{ o(t)+ = 6(0)6(0)] o(t) + s(H)6(0) >e€ 
u(t) = 
- X (oi +o wjeco)| -E< o(t) + 6(t)|6(1)} < € (3-34) 
E 2N 2N 


In the next chapter, several simulations will be presented in order 
to illustrate the effectiveness of the Bang-Bang control algorithms as compared 
to Proportional Navigation, in both the boost phase and in the terminal phase of 
the missile flight. 

3. Missile Dynamics 
The signal flow graphs for the missile dynamics are shown in Figure 


i) 
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Figure 15. Missile Dynamics 


In state-space, these signal flow graphs can be represented as 


Ke ee Ox. 0 0 

xe ORONO Oilex 1 O|] ucosy 

. |= 5 fc (3-35) 
Men Ca Oe | ye, O Oj} -usiny 

Va 00 0 Olly, 0 1 
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IV.VERTICAL LAUNCH MISSILE BOOST PHASE 
SIMULATION 


A. OVERVIEW 

In simulating the boost phase of the vertically launched missile, several 
assumptions can be made to reduce the complexity of the simulation to a level 
that won't detract from the concepts being presented: 

1. The missile is limited to Mach 4. 

2. The missile is limited to 30 g's of acceleration in the transverse and 
tangential directions combined. 

3. The speed of sound in air is constant at 340 m/s. 


4. The density of air is constant at 1 Kg/m3. 





5. The missile mass remains constant. 
6. The lift and drag forces experienced by the missile can be modeled by | 


equations (3-6) through (3-9). | 


B. MISSILE PARAMETERS | 
At launch, the following parameters define the state of the missile: 


The initial missile parameters are: 


Ly=4.2 m missile length 

M=225 Kg . missile mass 

L,pg=2.1 m length from center of gravity to tail 
Lpg=0.15 m length from center of gravity to center 


of pressure 


S=0.99 m2 missile reference lift surface 
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T=66800 N missile thrust 

I,,=2000 kg m2 missile moment of inertia 
o=1. kg/m3 density of air 

V max= 1460 m/s max missile velocity 


Additionally, the launch point will correspond to inertial reference: 
x(0)=0 m y(0)=0 m 
V,(0)=0 m/s V,(0)=0 m/s 


C. CANNISTER EGRESS 
It is not possible, nor is it desirable, to have the missile commence its tip- 
Over maneuver until it has egressed from the launch cannister and achieved 
sufficient distance from the launch platform to minimize hazards to personnel. 
Thus, the TVC actuators will have a zero degree deflection for the first 0.7 
seconds after launch. This will ensure sufficient missile altitude before tip-over 
begins. At time=0.7 seconds the missile position and velocity are: 
x(0.7)=0 m y(0.7)=71.8 m 
V,(0.7)=0 m/s V,(0.7)=204.5 m/s 
D. BOOST SIMULATION AND RESULTS 
The simulation was run with four values of the maximum thrust vector 
control (TVC) angle: 30°, 45°, 60°, and 90°. This was done to illustrate what 
effect a greater TVC angle will have on the missile trajectory and velocity. 


Figure 16 shows the missile trajectories for the four cases. 
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Figure 16. Missile Boost Phase Trajectories 


Analysis of the trajectories shows that as the maximum TVC angle 1s 
increased, the missile will achieve horizontal flight at lower altitudes than 
smaller TVC angles. However, at angles greater than 45°, there is a significant 
reduction in the downrange distance. This may be crucial if the missile is to be 
fired against a target that is close to the ship at launch. At TVC angles less 
than 30°, the thrust is directed toward achieving missile altitude to a larger 
extent than having the missile achieve horizontal flight in minimum time. This 
may be advantageous if the missile if fired against a high-altitude target, but not 
against a sea-skimming missile. Another consequence of having a large TVC 


angle is that much of the thrust is spent rotating the missile to horizontal at the 
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expense of missile kinetic energy, which results in a longer time for the missile 


to reach maximum velocity. This is shown graphically in Figure 17. 
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Figure 17. Missile Boost Phase Velocity Magnitude Profiles 


At the two highest TVC angles (60° and 90°), the first several seconds are 
spent turning the missile over while not attaining much forward velocity. As a 
result of a higher turning rate, higher angles of attack are generated, as shown 
in Figure 18. This is not desirable because at higher angles of attack, the drag 
forces that act to slow the missile down increase. Thus, it is advantageous 
from an aerodynamic point of view to keep the angle of attack as small as 
possible in order to minimize the drag force acting on the missile. 

Figure 19 shows the missile pitch angle as a function of time. It is of 
interest to note that there is not a great difference between total turning times 


between the 30° TVC angle and the 90° TVC angle. This is due to the counter- 
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acting lift and drag forces encountered by the missile, during high TVC angles, 
that are caused by high angles of attack. This indicates that increasing the 
maximum TVC angle beyond 45° may achieve a somewhat faster turning time, 
but at the expense of a great loss of kinetic energy and therefore velocity. This 
is of peak importance, for without sufficient velocity, a missile has little or no 


chance of successfully intercepting an incoming target, particularly a 


maneuvering target. 
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Figure 18. Missile Boost Phase Angle of Attack Profiles 
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Figure 19. Missile Boost Phase Pitch Angle Profiles 
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V. TERMINAL PHASE SIMULATIONS 


In order to compare the effectiveness of minimum-time ‘Bang-Bang' control 
as compared to proportional navigation control, several missile-target 
geometries were simulated. Each geometry was run with four control system 
configurations: 

1. proportional navigation using analytical values of line-of-sight rate and 
acceleration, which shall be referred to as ‘analytical proportional navigation’; 

2. proportional navigation using estimates of line-of-sight rate and 
acceleration, which shall be referred to as ‘estimated proportional navigation’; 

3. Bang-Bang control using analytical values of line-of-sight rate and 
acceleration, which shall be referred to as ‘analytical Bang-Bang;; 

4. Bang-Bang control using estimates of line-of-sight rate and acceleration, 


which shall be referred to as ‘estimated Bang-Bang; 


A. CASE ONE 
In the first set of simulations, the target will be headed directly toward the 


ship, flying a straight course with no maneuvers. The initial conditions are: 


Missile: 
x=0: z=1000 m; 
Vx=680 m/s (Mach 2) Vz=0; 
Target: 
x=7000 m; z=30 m 
Vx=-850 m/s (Mach 2.5); Vz=0; 
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Figure 20 shows the missile-target geometries using true line-of-sight rate 
and acceleration. Both proportional navigation and Bang-Bang control was able 
to hit the target. Note that the trajectory using Bang-Bang control flew a 
straighter course to the intercept point than the trajectory using proportional 
navigation control. This indicates that Bang-Bang control is able to guide the 
missile onto an intercept trajectory faster than the proportional-navigation 
controller, which is continuously changing the missile flight path in order to 
intercept the target. 

Figures 21 and 22 show control force as a function of time for the 
proportional navigation controller and the Bang-Bang controller, respectively. 


The proportional navigation controller required continuous control in order to 
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Figure 20. Case One: Missile/Target trajectories using true line 
of sight rate 
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correct the missile trajectory enough to enable the missile to intercept the 
target. It does not leave much reserve control to counter a maneuver away 
from the missile by the target. Conversely, the Bang-Bang controller applied 
maximum control effort until the line of sight rate became zero, then shut off. 
This allows the missile greater flexibility in case the target maneuvers in any 
direction. If the linear zone were not incorporated into the Bang-Bang 
controller, however, the controller would have chattered between maximum 
positive and negative control force instead of shutting off. This would have 
decreased the missile velocity unnecessarily just when the missile needs 


maximum kinetic energy for maneuverability. 
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Figure 21. Case One: Control force versus time using proportional 
navigation control with true line of sight rate 


Figure 23 shows the trajectories using estimates of line of sight rate and 
acceleration. It is apparent here that the proportional navigation controller was 
not able to react fast enough in order to hit the target. This is due to the delay 
encountered while the observer filter was not matched up with the analytical 
values, as shown in Figure 24. By the time the estimated value of line of sight 
rate and acceleration matched the analytical values, the target was at a point 


where proportional navigation could not generate enough of a control signal fast 
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enough to enable target intercept. This illustrates the importance of being able 
to detect, virtually instantaneously, changes in target velocity and direction. The 
Bang-Bang controller also experienced the same delay before the proper 
control was applied, as shown in Figure 25, but since it immediately applied 


maximum control, the missile was still able to intercept the target. 
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Figure 22. Case One: Control force versus time using Bang-Bang 
control with true line of sight rate 
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Figure 23. Case One: Missile/Target trajectories using estimated 
values for line of sight rate 
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Figure 24. Case One: Control force versus time using proportional 
navigation control with estimated values for Ime of sight rate 


acceleration (m/s 4 


lS 2 
time (sec) 





Figure 25. Case One: Control force versus time using Bang-Bang 
control with estimated values for line of sight rate 


Figure 26 compares the trajectories of Bang-Bang control using analytical 
and estimated values of line of sight rate and acceleration. Because the 
estimated values took approximately 0.25 seconds to match the analytical 
values, there is a period at the start of the simulation where the estimated 
trajectory has positive control applied to it, while the analytic trajectory has 
negative control applied. After the analytic and estimated values matched up, 
they both were able to drive the line of sight rate to zero in plenty of time to 


intercept the target. 
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Figure 26. Case One: Missile/Target trajectories using Bang-Bang 
control, comparing true and estimated values for line of sight rate 


B. CASE TWO 
In the second set of simulations the target will be headed directly toward 


the ship, but will initiate a 12 g vertical maneuver at problem start. The initial 


conditions are: 


Missile: 
X=0; z=1000 m 
Vx=680 m/s (Mach 2) Vz=0; 
Target: 
x=7000 m z=30 m 
Vx=-850 m/s (Mach 2.5); Vz=0) 
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Figure 27 shows the simulation results of the analytical proportional 
navigation and analytical Bang-Bang controllers. The Bang-Bang controlled 
missile was able to adjust to the 12 g vertical maneuver of the target in order to 
achieve intercept. The proportional navigation controlled missile was not able 
to keep up with the high g maneuver of the target, and missed the target by 
about 80 meters. 

Figure 28 shows the control force as a function of time for the analytical 
proportional navigation controlled missile. Note that the controller did not 
saturate until well into the simulation. This indicates that if more control had 
been applied sooner, the missile could have intercepted the target. This is 


illustrated by Figure 29, which shows the control of the analytical Bang-Bang 
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Figure 27. Case Two: Missile/Target trajectories using true line of 
sight rate 
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Figure 28. Case Two: Control force versus time using Proportional 
Navigation control with true line of sight rate 


= 


acceleration (m/s‘2) 


time (sec) 





Figure 29. Case Two: Control force versus time using Bang-Bang 
control with true line of sight rate 


controlled missile. Neither controller went to zero as the Bang-Bang controllers 
of Case One did. This is because the target is maneuvering, constantly 
changing the line of sight angle. 

The simulations using estimated values for line of sight rate and 
acceleration are shown in Figure 30. As was the case with the analytical 
expressions, the Bang-Bang controller using estimates was still able to achieve 
intercept, while the proportional navigation controller missed the target. Note 


that the miss distance for the proportional navigation controller in Figure 30 is 
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less than the miss distance achieved with the analytical expressions. This is 
because the target turned toward the missile, and the delay caused by the 
observer prevented the missile from beginning to dive after the target as 
quickly, thus leaving less of an altitude difference when control is applied. If the 
target were to turn away from the missile, the estimates of line of sight rate and 
acceleration would be worse than those obtained using analytical values. 
Figures 31 and 32 are the control versus time graphs for the estimated 
proportional navigation controller and the estimated Bang-Bang controller. The 


delay encountered by both of the controllers is about the same as those 


encountered in Case One. 
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Figure 30. Case Two: Missile/Target trajectories using estimated 
values for line of sight rate 
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Figure 31. Case Two: Control force versus time using proportional 
navigation control with estimated values for line of sight rate 
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Figure 32. Case Two: Control force versus time using Bang-Bang 
control with estimated values for line of sight rate 
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VI. CONCLUSIONS AND RECOMMENDATIONS 


This work has shown that the minimum-time (Bang-Bang) controller is an 
effective algorithm for missile control, in both the boost and terminal phase of 
tactical missile flight. It is particularly effective when the target has a speed 
advantage over the missile or when the target is maneuvering. 

As is the case with virtually all tactical missiles, the accuracy of the target's 
measured position, velocity, and acceleration vectors are paramount in 
accurately predicting the parameters with which the missile is controlled with a 
Bang-Bang minimum time controller. 

Areas for future study will include: 

1. Simulating the entire missile flight from launch through intercept. The 
focus of these studies will be to analyze the effects of missile velocity, tip-over 
altitude, and target speed advantage on missile performance. 

2. Developing more complex models in order to better understand what 
the effects of a Bang-Bang controller will have on a particular missile system. 
Such studies will use empirical aerodynamic and physical data for a particular 
missile, such as SM-II Block IV or Vertically-Launched Sea Sparrow, and will 
entail developing a 5 or 6 DOF computer model. 

3. Investigating the effects of noise on measurement of line of sight and 
its time derivatives, to the extent of the effect noise has on a Bang-Bang 
controller. 

4. Investigating the effects of plant modelling errors on system sensitivity 


when using Bang-Bang controllers. 
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APPENDIX 1- BOOST PHASE MATLAB PROGRAM 


% VLS BOOST PHASE MISSILE SIMULATION 
% T B MULL 

Jo REVISED 22 APRIL 1992 

clear;clg; 

Jo Initial Missile Parameters 

Delta=[30 45 60 90}; 

for j=1:length(Delta) 


L_m=3.6; Yo L_m= missile length in meters 

M=227; Yo M= missile mass in kilograms 

L_gt=1.8; %. L_gt= length from center of gravity to tail 

Peep=0215. % WL_gp= length from center of gravity to center 
Jo of pressure 

S=0.99; % S 1s lift surface in square meters 

T=66800; % T= missile thrust in newtons 

Icg=2000; %o. Icg is the missile moment of inertia (kg m‘’2) 

Tho=4R00; %o rho is the density of air in kg/m‘3 

Vmax=340*4; Yo Vmax is the max missile velocity in m/s 

di= 02, % dt is simulation time step size 


delta_max=Delta(j); 


Op averse oe ee 
Jo Initial Conditions 
tfinal=10.00; % Final simulation time. 
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kmax=tfinal/dt; 
time=0.0; 
x=0;xdot=0;xddot=0; 
y=0;ydot=0;yddot=0; 


gamma=90.0; 
gammadot=0; 

theta=90.0; 

thetadot=0; 
X_theta=[thetadot;theta]; 
X_x=[xdot;x]; 
X_y=lydot;y]; 

A=[0 0;1 0]; 


B=[1;0]; 
[phi,del]=c2d(A,B,dt); 
alpha=theta-gamma; 


dtr=pi/180; 


% kmax is number of time increments 

Zo initial time 

% initial missile states (in x) 

% initial missile states (in y) 

J (these initial y states correspond to missile 
Yo position and velocity after clearing the 

% cannister) 

Yo gamma is inertial to velocity vector 

% in degrees (will convert to radians in loop) 


% theta is inertial to missile head 


% A and B matrices for determining velocity and 
% position of x,y, and theta base on acceleration 


% commands 


% converts A,B to discrete for simulation 
% alpha is Angle-of-attack 


% converts degrees to radians 


V_m(1)=sqrt(X_x(1,1)42+X_y(1,1)*2); % V_m is missile velocity magnitude 


F_drag=0; 
F_lift=0; 
Cd=0; 
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for i=1:kmax-1 


Jo****aloorithm for determining rudder deflection angle*************##e** 
9p BA CEE CIC CIC ICICI ACI ICICI ICI AIC A ACCA AC ACCA SEO A CoA ER ICA HCA EA EA 
% test for cannister egress(time<0.6) 
if time(i)<0.7 

delta(i)=0; 
else 

Theta=X_theta(2,1); 

Thetadot=X_theta(1,1); 
N=T*L_gt/Icg*sin(delta_max*dtr); 
test(i)=Theta+(Thetadot)*abs(Thetadot)/(2*N); 

eps=0.5; 


if abs(test(i))>eps 
delta(i)=delta_max*sign(test(i)); 
else 
delta(i)=delta_max/eps*test(i); 
end; 
end; 
Op 707 Ae A oe A oe Ae oe 2 ae Ae 2 oe 6 oe oe A ae A ae 2 G2 Ae OE Ae 2 AC OR 2K EO oR oe oe eo EO eo A eB oe oe oe oA aS oh oa oe 
F_lift(:,1)=0.5*rho*V_m(1i)*2*S*0. 1 *alpha(i); 
Cddi)=(2*T/(M *rho* Vmax%2*S)); 
F_drag(:,1)=Cd(1)*0.5*rho* V_m(i)*2*S; 
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%F_drag(: ,1)=0; 
xddot(i)=-F_lift(i)/M*sin(gamma(i)*dtr)+T/M*cos(gamma(i)*dtr+delta(i) *dtr)... 
-F_drag(:,i)*cos(gamma(i)*dtr); 
yddot(i)=F_lift(1)/M*cos(gamma(i)*dtr)+T/M*sin(gamma(i)*dtr+... 
delta(i)*dtr)-F_drag(:,1)*sin(gamma(i)*dtr); 
thetaddot(1)=-F_lift(i)*L_gp/Icg*cos(alpha(i)*dtr)-T*L_gt* sin(delta(i)*dtr)/Icg; 


X_theta(:,i+1)=phi*X_theta(:,1)+del*thetaddot(i); 
X_x(:,1+1)=phi*X_x(:,1)+del*xddot(i); 
X_y(:,1+1)=phi*X_y(-,1)+del*yddot(1); 


gamma(i+1)=atan2(X_y(1,i+1),X_x(1,i+1))/dtr; 
% X_y=[ydot;y]; 
alpha(i+1)=X_theta(2,i+1)-gamma(i+1); 
delta(i+1 )=delta(i); 
V_m(i+1)=sqrt(X_x(1,1+1)42+X_y(1,i+1)*2); 
time(i+1)=time(i)+dt; 

end; J END OF LOOP 

J AAA OCC CACAO CCCI ASAIO IA IEA IA A AGK 
eval(('"X_x' num2str(j),'=X_x;']); 
eval(("X_y',num2str(j),'=X_y;']); 
eval({'V_m’',num2str(j),,=V_m;']); 

eval([‘alpha' ,num2str(j),'=alpha;']); 
eval(('X_theta' num2str(j),'=X_theta;']); 


delta_max 


> 


vm=V_m(it1) 
X=X_x(2,i+1) 
Y=X_y(2,i+1) 
Alpha=alpha(i+1) 
Theta=X_theta(2,1+1) 
end; 

Boost] plot2a 


save bocstdata 


clg; 
z=length(X_x(2,:)); 
for k=1:length(Delta) 
eval({'X_x=X_x',num2str(k),’';']); 
eval({"X_y=X_y',num2str(k),’;']); 
axis({[O 9000 0 2000)); 
plot(X_x(2,:),X_y(2,:));grid; 
hold on 
if k==1 
plot(X_x(2,1:2/dt:z),X_y(2,1:2/dt:z),'o’) 
plot(X_x(2,z),X_y(2,z),'0') 
elseif k== 
plot(X_x(2,1:2/dt:z),X_y(2,1:2/dt:z),'o') 
plot(X_x(2,z),X_y(2,z),'0’) 
else 


plot(X_x(2,1:2/dt:z),X_y(2,1:2/dt:z),'o') 


a 





plot(X_x(2,z),X_y(2,z),'0') 
end; 
end; 
hold off 
title( Missile Boost Phase- Bang-Bang control’); 
xlabel( Horizontal Range (m)’),ylabel(Altitude (m)’); 
%text(0.6,0.8,_ 30° TVC angle’,'sc’) 
pause;clg; 
axis( normal’) 
%o 
for k=1:length(Delta) 
eval(("V_m=V_m' ,num2str(k),';']); 
plot(time, V_m); 
hold on; 
end; 
hold off 
title("V_m versus time’);grid; 
xlabel(‘Time (sec)’),ylabel('Velocity (m/s)’); 
%o 
pause;clg; 
for k=length(Delta):-1:1 
eval([‘alpha=alpha' num2str(k),';']}); 
plot(time,alpha),grid; 
hold on; 


end; 
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hold off 

title(Plot of Alpha versus time’) 

xlabel(‘Time (sec)'); 

ylabel(’'Angle of Attack (degrees)’); 

pause;clg; 

Jo 

axis([O 5 -20 100]); 

for k=1:length(Delta) 
eval(['X_theta=X_theta',num2str(k),';']); 
plot(time,X_theta(2,:)), grid; 
hold on; 

end; 

hold off 

title('Plot of Theta versus time’) 

xlabel(‘Time (sec)'); 

ylabel(‘Theta (degrees)'); 


axis(‘normal’) 
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APPENDIX 2- TERMINAL PHASE MATLAB PROGRAM 


% terminal.m 


Jo TERMINAL MISSILE-TARGET ENGAGEMENT 
Jo TIM MULL 
% 24 APRIL 1992 1130 


% This model does not take drag or lift into consideration 

% -uses two different methods to compute sigma derivatives: 

% 1. calculates sig_dot, sig_ddot using analytical expressions 

% 2. uses an observer to estimate sighat,sigdothat,sigddothat 

% 3. controls using sigmadothat,sigmaddothat 

% A. utilized bang bang on sigmadot,sigmaddot with a small linear region 
% in order to alleviate controller chatter once lock has been achieved. 


clear;clg; 


cle; 

disp([’ I); 

for k=1:4 

disp(['Control= ',num2str(k)]); 
control=k; 


dt=0.01;, 


Op BEER AK KKK K Initial Conditions FE OG OE 2g 2 RR OG OIE Og 2 2 RS OE 2 oe OE 2 OK oo OK OK OK OOK OK OK OR OK OK KK KK 


a) 


xm=0; 
vxm=2*340; 
axm=0; 
ym=1000; 
vym=0; 
aym=0; 


xt=5000; 
vxt=-2.5*340; 
axt=0; 

yt=30; 

vyt=0; 

ayt=0; 


time=-0; 


%observer stuff 


pl=-100; 
p2=-100; 
p3=-100; 


% missile horizontal position in meters 
% missile horizontal velocity in m/s 
% missile horizontal acceleration in m/s%2 
% missile vertical position in meters 
% missile vertical velocity in m/s 


% missile vertical acceleration in m/s“2 


% target horizontal position in meters — 
% target horizontal velocity in m/s 
% target horizontal acceleration in m/s*2 
% target vertical position in meters 
% target vertical velocity in m/s 


% target vertical acceleration in m/s%2 


K=conv({1 -pl],conv({1 -p2],[1 -p3])); 


k1=K(4); 
k2=K(3); 
k3=K(2); 
X=[0;0;0]; 
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for i=1:1200 Zoloop | 
%if rem(i,100)==0 
Joclg; 
Josubplot(221),plot(xm,ym) 
Josubplot(222),plot(time(1 :length(U)),U) 
% end; 
gamma(i)=atan2(vym(i),vxm(i)); % gamma is angle of missile velocity 
% vector with respect to inertial ref 


Vm(i)=sqrt(vxm(i)*2+vym(i)*2); % Vm is missile velocity magnitude 


A(i)=xt(1)-xm(i); % X is relative horizontal position 
Y(i)=yt(i)-ym(1); % Y is relative vertical position 
VX(i)=vxt(i)-vxm(i); % VX is relative horizontal velocity 
VY(@)=vyt(i)-vym(i); % VY is relative vertical velocity 
AX(1)=axt(i)-axm(1); % AX is relative horizontal acceleration 
AY(i)=ayt(i)-aym(i); % AY is relative vertical acceleration 


RG)=sqrt(X(i)*2+ Y(1)*2); 
Rdot(ij=sqrt(VX(1)*2+V Y (i)*2); 
% sigma,sigmadot,sigmaddot are the line 
% of sight between the missile and the 
% target, and its Ist 2 time derivatives 
Of 8 AAAI NOTE, 6 Ca Eo ok ACE ECE A EA EA A EE EC EI 
% This simulation is using analytical expressions for sigma * 
%o sigmadot, and sigmaddot. Since the target's velocity and * 


% acceleration cannot be measures directly, a more realistic * 
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% approach is to use an observer to estimate sigmadothat and * 

%. sigmaddothat. * 

Of 3 5 CACC ACR CR ACA CR ACR ACE CA CA CE AAC AC ACE fC AC A A A A 
sigma(i)=atan2(Y(1),X(i)); 
sigmadot(i)=(X(i)*V Y (i)- Y 1) * VX(1))/RG)2; 
sigmaddot(i)=(X(i)*A Y(i)- YG) *AX())/RG)2... 

-2*(X()* VX()+Y(1)* VY (1) *(X)* V ¥G)- YQ)*V.XG))/RG)"4; 
OEE EAC CCC CACC ACI ACACIA AC ACACA ICI CAC A A CCA AC EA HC 
Yo 3rd order observer to obtain sigma.dot.hat,sigma.ddot.hat 

A=[-k3 1 0 
-k2 -2*Rdot(i)/R@) 1 
-k1 O O}; 
B=[ k3 
k2 
kL]; 
[phi,del]=c2d(A,B,dt); 


X(:,1+1)=phi*X(:,1)+del*sigma(1); 
sigmadothat(i)=X(2,i); 
sigmaddothat(i)=X(3,i); 


A FE HE AG AG A He Ae 2 ie Ae 2 ee ee 6 2 ie He 2 fe Ae 2 Ae He 2 Ae Ae a Ae oe Ae oR ee Oe RE EE ee Ae ee Ee oo 


Yo CONTROL SECTION 
Yo 
Yo PROP NAV CONTROL 
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if control== % if loop 1 
U@=Vm(i)*4*sigmadot(i); 
if abs(U(i))>15*9.81 % if loop 2 
U(@)=15*9.81*sign(U(i)); 
end; % end if loop 2 
Jo 
Jo PROP NAV CONTROL based on sigmadothat 
elseif control== % if loop 1 
U@=Vm()*4*sigmadothat(i); 
if abs(U(1))>15*9.81 % if loop 2 
U()=15*9.81*sign(U(i)); 
end; % end if loop 2 
Jo 
elseif control== % if loop 1 
Jo Bang-Bang control based on sigmadothat, sigmaddothat 
N=15*9.81; % N is maximum normal acceration in Newtons 
test(i)=sigmadothat(1)+(sigmaddothat(i)*abs(sigmaddothat(i)))/(2*N); 
eps=0.001; 
if abs(test(i))>eps % if loop 3 
U@=N*sign(test()); 
else 
U(G)=N/eps*test(i); 
end; % end if loop 3 
Jo 


elseif control==4 
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To Bang-Bang control based on sigmadot, sigmaddot 
N=1579°S1- % N is maximum normal acceration in Newtons 
test(1)=sigmadot(i)+(sigmaddot(i)*abs(sigmaddot(1)))(2*N); 
eps=0.001; 
if abs(test(i))>eps 

U(@=N*sign(test(1)); 
else 
U(G)=N/eps*test(i); 
end; 
end; 
Jo 
if rem(i,50)==0 
disp('working’); 
end; 


% ee Mie missile States ®#EEEE EE EER KEE EKA EKA AL A A Ee 


% Let Xmsl=[xm(i) State space representation of missile states 
% vxm(i) 

%o ym(i) 

To vym(i)]; 


Xmsl(:,1)=[xm(i) 
vxm(1) 
ym(1) 
vym(1)]; 


Am=[0 | 00 
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0000 
0001 
000 0]; 


Bm=[0;-sin(sigma(i));0;cos(sigma(i))]; % control force is perp to sigma!! 
[phim,delm]=c2d(Am,Bm,dt); % convert continuous state-space to discrete 


% update missile states 
Xmsl(:,1+1)=phim*Xms\l(:,1)+delm*U(1); 


xm(i+1)=Xmsl(1,i+1); 
vam(i+1)=Xms|(2,i+1); 
axm(i+1)=-sin(sigma(i))*U(); 
ym(i+1)=Xmsl(3,i+1); 
vym(it+1)=Xmsl(4,i+1); 
aym(i+1)=cos(sigma(i))*U(i); 


% jae **ypdate target States 6 ee Eee TEE EE EEE ee eee 


% Let Xtet=[xt(i) 


Jo vxt(1) 

Jo yt(1) 

7 vyt@)]; 

Atgt(:,=[xt(i) 
vxt(1) 
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yt(i) 
vyt(i)]; 


At=[0 100 

0000 

0001 

O10 0-0], 
Bt=[0;axt(i);0;ayt(i)]; 


[phit,delt]=c2d(At,Bt,dt); 
Xtegt(:,i+ 1)=phit*Xtgt(:,1)+delt*1; 


xtit+1)=Xtgt(1,i+1); 
vxt(it+1 )=Xtet(2,i+1); 
axt(it+1)=axt(i); 
ytG+1)=Xtgt(3,i+1); 
vytGitl)=Xtgt(4,1+1); 
ayt(it+1)=ayt(1); 
time(it+ | )=time(i)+dt; 


Totes eae Check (Of Cra 


X(i+1 )=xt(it+1)-xm(i+1); % X is relative horizontal position 


Y Gi+1)=ytG+1)-ym(i+1); % Y is relative vertical position 
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R(i+1)=sqrt(X(it+1)*2+ Y(i+1)%2); 


Jo**KAKEE subroutine to interpolate the 2 closest ranges to find CPA 
if RG+1)>RQG) 
xms=[ym(i-1:1+1);vym(i-1:1+1);xm(i-1:1+1);vxmG-1:1+1)]; 
xts=[ytG-1:1+1);vytQi-1:1+1);xtQ-1:1+1);vxtG-1:1+1)]; 
rl = interp(xms(:,1:2),xts(:,1:2)); 
r2 = interp(xms(:,2:3),xts(:,2:3)); 
Rmin=min((rl r2]); 
break; 
end; 
end; 
j=k; 
eval({'xm’',num2str(j),’ =xm;']); 
eval(['ym',num2str(j),' =ym;']); 
eval(('xt',.num2str(j),’ =xt;']); 
eval([‘yt ,num2str(j),’ =yt;']); 
eval({'Rmin',num2str(j),, =Rmin;'}); 
eval({'U',num2str(j),' =U;']); 
eval(['save data',num2str(j)]); 


end; 
clg; 


clear 


load datal 
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save dataS xml yml xtl ytl 


clear 

load data2 

xm2=xm; 

ym2=ym; 

XIZ=Xt- 

yt2=yt; 

save data5 xml ym] xtl ytl xm2 ym2 xt2 yt2 


clear 

load data3 

xm3=xm; 

ym3=ym; 

Kio=KG 

yt3=yt; 

save dataS xm] ym] xtl ytl xm2 ym2 xt2 yt2 xm3 ym3 xt3 yB 


clear 

load data4 
xm4=xm; 
ym4=ym; 
X4=xt 


yt4=yt; 
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save data5 xml ym! xtl ytl xm2 ym2 xt2 yt2 xm3 ym3 xt3 yt3 xm4 ym4 xt4 
yt4 


clear 


% This routine takes data computed in terminal.m and samples it every 


% 8 data point in order to make the graphs less memory intensive 


clear 


load data5 


xmla=[];ym1a=[];xtla=[];ytla=[];xm2a=[]; ym2a=[];xm3a=[];ym3a=[]; 
xm4a=[];ym4a=[]; 

k=min({length(xm1 ),length(xm2),length(xm3),length(xm4)]); 

for i=1:k 


if rem(i,8)==0 


xmla=[xmla xml(i)]; 
ymla=[ymla yml(i)]; 
xtla=[xtla xtlq)]; 
ytla=[ytla ytl@)]; 
xm2a=[xm2a xm2(i)]; 
ym2a=[ym2a ym2(i)]; 
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xm3a=[xm3a xm3(i)]; 

ym3a=[ym3a ym3()J; 

xm4a=[xm4a xm4(i)]; 

ym4a=[ym4a ym4(i)]; 
end; 
end; 
axis({O 5000 0 1200]); 
plot(xm1la,ym1a,'-w',xtla,ytla,'-w',xm4a,ym4a,'-w ) 
pause; 


plot(xm2a,ym2a,'-w',xtla,ytla,'-w',xm3a,ym3a, -w') 
cleg;clear 


load data] 
ul=U- 

load data2 
uZ=U;: 

load data3 
u3=U; 

load data4 
u4=U; 

save control ul u2 u3 u4 
clear 

load control 


dt=0.01; 
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n=min((length(u1),length(u2),length(u3),length(u4)]); 
tfinal=(n-1)*dt; 

time=0.00:dt:tfinal; 
subplot(211),plot(time,u1(1:n),'-w'),titleCcontrol=1’); 
subplot(212),plot(time,u2(1:n),'-w’),title(control=2’); 
pause;clg; 
subplot(211),plot(time,u3(1:n),'-w'),titleCcontrol=3’'); 
subplot(212),plot(time,u4(1:n),'-w’),title(control=4'); 
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